Stevens East

2013

plot.gps.grids(corn.2013.dat, swath.width = 12.6, plots=0,residuals=TRUE) 

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

2014

plot.gps.grids(soybean.2014.dat,plots=0,residuals=TRUE) 

2015

plot.gps.grids(corn.2015.dat,swath.width = 12.6, plots=0,residuals=TRUE) 

2016

plot.gps.grids(soybean.2016.dat,plots=0,residuals=TRUE) 

2017

plot.gps.grids(corn.2017.dat,swath.width = 12.6,plots=0,residuals=TRUE) 

2018

plot.gps.grids(soybean.2018.dat,plots=0,residuals=TRUE) 

2020

plot.gps.grids(corn.2020.dat,swath.width = 12.6,plots=0,residuals=TRUE) 

2021

plot.gps.grids(soybean.2021.dat,plots=0,residuals=TRUE) 

Variograms

Yield Monitor

Full Field

corn.2013.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=corn.2013.dat)
plot(corn.2013.var)

soybean.2014.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=soybean.2014.dat)
plot(soybean.2014.var)

corn.2015.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=corn.2015.dat)
plot(corn.2015.var)

soybean.2016.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=soybean.2016.dat)
plot(soybean.2016.var)

corn.2017.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=corn.2017.dat)
plot(corn.2017.var)

soybean.2018.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=soybean.2018.dat)
plot(soybean.2018.var)

corn.2020.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=corn.2020.dat)
plot(corn.2020.var)

soybean.2021.var <- variogram(Yield~1, 
                          locations=~Easting + Northing,
                          alpha=c(0,45,90),
                          data=soybean.2021.dat)
plot(soybean.2021.var)

Functional models for Random Fields

LOESS

#when using VRYIELDRY
#spans <- c(0.05, 0.1, 0.2, 0.4, 0.6, 0.8)

#note - at 2^(-6) the edges are problematic
spans <- 2^(-5:0)

LOESS <- vector(length(spans),mode='list')
Variograms <- vector(length(spans),mode='list')
Residuals <- vector(length(spans),mode='list')
Estimates <- vector(length(spans),mode='list')
loess.dat <- NULL
idx <- 0
for(span in spans) {
  idx <- idx + 1
  tmp.dat <- corn.2015.dat[,c('Northing','Easting','Yield')]
  tmp.loess <- loess(Yield ~ Easting + Northing, span=span ,data=tmp.dat)
  tmp.dat$Yield <- predict(tmp.loess)
  tmp.dat$residuals <- residuals(tmp.loess)
  tmp.dat$Span = span
  loess.dat <- rbind(loess.dat, tmp.dat)
  Estimates[[idx]] <- tmp.dat
  tmp.var <- variogram(Yield~1, 
                          locations=~Northing+Easting,
                          alpha=c(0,45,90),
                          data=tmp.dat)
  Variograms[[idx]] <- tmp.var
  #plot(tmp.dat$Yield)
  tmp.var <- variogram(residuals~1, 
                          locations=~Northing+Easting,
                          alpha=c(0,45,90),
                          data=tmp.dat)
  Residuals[[idx]] <- tmp.var
  #plot(tmp.var,main=paste('residuals',span))
}

#loess.dat <- rbind(loess.dat, tmp.dat)

ggplot(loess.dat, aes(Easting,Northing)) + 
geom_point(aes(colour = Yield),size=1) + 
scale_colour_gradient2(low=vermillion, mid=yellow, high=blue,midpoint = mean(loess.dat$Yield)) +
#scale_colour_gradient(low=cbPalette[8], high=cbPalette[5]) +
labs(colour = "LOESS", x="X (m)", y="Y (m)") + facet_wrap(~ Span)

ggplot(loess.dat, aes(Easting,Northing)) + 
geom_point(aes(colour = residuals),size=1) + 
scale_colour_gradient2(low=vermillion, mid=yellow, high=blue,midpoint = 0) +
#scale_colour_gradient(low=cbPalette[8], high=cbPalette[5]) +
labs(colour = "Residuals", x="X (m)", y="Y (m)") + facet_wrap(~ Span)

plot(Variograms[[1]],main=paste('LOESS',spans[1]))

plot(Variograms[[2]],main=paste('LOESS',spans[2]))

plot(Variograms[[3]],main=paste('LOESS',spans[3]))

plot(Variograms[[4]],main=paste('LOESS',spans[4]))

plot(Variograms[[5]],main=paste('LOESS',spans[5]))

plot(Variograms[[6]],main=paste('LOESS',spans[6]))

plot(Residuals[[1]],main=paste('Residuals',spans[1]))

plot(Residuals[[2]],main=paste('Residuals',spans[2]))

plot(Residuals[[3]],main=paste('Residuals',spans[3]))

plot(Residuals[[4]],main=paste('Residuals',spans[4]))

plot(Residuals[[5]],main=paste('Residuals',spans[5]))

plot(Residuals[[6]],main=paste('Residuals',spans[6]))

library(moments)
par(mfrow=c(2,3))
for(i in 1:6) {
  hist(Estimates[[i]]$Yield,
       main=paste('LOESS',spans[i]),
       xlab=paste("sk=",format(skewness(Estimates[[i]]$Yield),digits=4),
                  "ku = ",format(kurtosis(Estimates[[i]]$Yield),digits=4))
       )
}

library(moments)
par(mfrow=c(2,3))
for(i in 1:6) {
  hist(Estimates[[i]]$residuals,
       main=paste('Residuals',spans[i]),
       xlab=paste("sk=",format(skewness(Estimates[[i]]$residuals),digits=4),
                  "ku = ",format(kurtosis(Estimates[[i]]$residuals),digits=4))
       )
}

Borrowed from 7.grid_cells.Rmd (Workshop)

grid.dists <- as.matrix(dist(cbind(corn.2015.dat$Easting,corn.2015.dat$Northing)))
grid.dists <- 1/grid.dists
diag(grid.dists) <- 0
library(ape)

Moran.tbl <- data.frame(Source=rep(0,7),
                        Expected=rep(0,7),
                        MoranI.Y=rep(NA,7),
                        SD.Y=rep(NA,7),
                        p.Y=rep(NA,7),
                        MoranI.Res=rep(NA,7),
                        SD.Res=rep(NA,7),
                        p.Res=rep(NA,7))

yield.I <- Moran.I(corn.2015.dat$Yield, grid.dists)
Moran.tbl$Source[1] <- 0
Moran.tbl$MoranI.Y[1] <- yield.I$observed
Moran.tbl$Expected <- yield.I$expected
Moran.tbl$SD.Y[1] <- yield.I$sd
Moran.tbl$p.Y[1] <- yield.I$p.value
Moran.tbl$MoranI.Res[1] <- 0
Moran.tbl$p.Res[1] <- 1
for(i in 1:6) {
  yield.I <-Moran.I(Estimates[[i]]$Yield, grid.dists)
  res.I <-Moran.I(Estimates[[i]]$residuals, grid.dists)
  Moran.tbl$Source[i+1] <- spans[i]
  Moran.tbl$MoranI.Y[i+1] <- yield.I$observed
  Moran.tbl$SD.Y[i+1] <- yield.I$sd
  Moran.tbl$p.Y[i+1] <- yield.I$p.value
  Moran.tbl$MoranI.Res[i+1] <- res.I$observed
  Moran.tbl$SD.Res[i+1] <- res.I$sd
  Moran.tbl$p.Res[i+1] <- res.I$p.value
}
Moran.tbl

##    Source      Expected  MoranI.Y        SD.Y p.Y MoranI.Res      SD.Res p.Res
## 1 0.00000 -0.0002365744 0.1754390 0.001198324   0 0.00000000          NA     1
## 2 0.03125 -0.0002365744 0.2584453 0.001198295   0 0.01153639 0.001197848     0
## 3 0.06250 -0.0002365744 0.2778514 0.001198291   0 0.02583169 0.001197959     0
## 4 0.12500 -0.0002365744 0.2807251 0.001198274   0 0.04409356 0.001198034     0
## 5 0.25000 -0.0002365744 0.3059531 0.001198277   0 0.05066586 0.001198045     0
## 6 0.50000 -0.0002365744 0.3229100 0.001198294   0 0.05738791 0.001198084     0
## 7 1.00000 -0.0002365744 0.3536829 0.001198200   0 0.06810397 0.001198139     0

par(mfrow=c(1,3))
plot(MoranI.Y ~ Source,data=Moran.tbl)
plot(MoranI.Res ~ Source,data=Moran.tbl)
plot(MoranI.Res ~ MoranI.Y,data=Moran.tbl)

#fit<-lm(Yield ~ bs(Northing + Easting),data = corn.2015.dat)

Yield1.lm <- lm(Yield ~ poly(Northing, Easting, degree=5), data=corn.2015.dat)

The response surface library rsm has some useful plots to help us visualize the trend surface

library(rsm)
par(mfrow=c(1,3))
image(Yield1.lm, Northing ~ Easting)
contour(Yield1.lm,Northing ~ Easting, image = TRUE)
persp(Yield1.lm, Northing ~ Easting, zlab = "Yield, Poly 5")

#image(x=loess.dat$Northing, y=loess.dat$Easting, z=loess.dat$Yield)
#contour(x=loess.dat$Northing, y=loess.dat$Easting, z=loess.dat$Yield, image = TRUE)

B-splines

Perperoglou (2019) note that the bs spline function in R may produce spurious results when used with 2-dimensional data.

Natural Splines

Generalized Additive Models

gamlss

mgcv

Adapted from https://stackoverflow.com/questions/16554690/interpolate-in-a-3-dimensional-spline-in-r

library(mgcv)

## Loading required package: nlme

## This is mgcv 1.8-41. For overview type 'help("mgcv-package")'.

b1 <- gam(Yield ~ s(Easting,Northing),data=corn.2015.dat)
vis.gam(b1)
title("t.p.r.s")

corn.2015.dat$Spline <- predict(b1)

par(mfrow=c(1,3))
image(b1, Northing ~ Easting)
contour(b1,Northing ~ Easting, image = TRUE)
persp(b1, Northing ~ Easting, zlab = "Yield, Spline")

ggplot(corn.2015.dat, aes(Easting,Northing)) + 
geom_point(aes(colour = Spline),size=1) + 
scale_colour_gradient2(low=vermillion, mid=yellow, high=blue,midpoint = mean(corn.2015.dat$Spline)) +
#scale_colour_gradient(low=cbPalette[8], high=cbPalette[5]) +
labs(colour = "Residuals", x="X (m)", y="Y (m)")

model 2 with te() smooths

b2 <- gam(Yield ~ te(Easting,Northing),data=corn.2015.dat)
vis.gam(b2)
title("tensor product")

par(mfrow=c(1,3))
image(b2, Northing ~ Easting)
contour(b2,Northing ~ Easting, image = TRUE)
persp(b2, Northing ~ Easting, zlab = "Yield, Poly 5")

model 3 te() smooths specifying margin bases

b3 <- gam(Yield ~ te(Easting,Northing,bs=c("tp", "tp")),data=corn.2015.dat)
vis.gam(b3)
title("tensor product")

par(mfrow=c(1,3))
image(b3, Northing ~ Easting)
contour(b3,Northing ~ Easting, image = TRUE)
persp(b3, Northing ~ Easting, zlab = "Yield, TE")

par(mfrow=c(1,3))
contour(b1,Northing ~ Easting, image = TRUE)
contour(b2,Northing ~ Easting, image = TRUE)
contour(b3,Northing ~ Easting, image = TRUE)

gam1 <- gam(Yield ~ s(Easting,Northing),data=corn.2013.dat)
gam2 <- gam(Yield ~ s(Easting,Northing),data=corn.2015.dat)
gam3 <- gam(Yield ~ s(Easting,Northing),data=soybean.2016.dat)
gam4 <- gam(Yield ~ s(Easting,Northing),data=corn.2017.dat)
gam5 <- gam(Yield ~ s(Easting,Northing),data=soybean.2018.dat)

par(mfrow=c(2,3))
contour(gam1,Northing ~ Easting, image = TRUE)
contour(gam2,Northing ~ Easting, image = TRUE)
contour(gam3,Northing ~ Easting, image = TRUE)
contour(gam4,Northing ~ Easting, image = TRUE)
contour(gam5,Northing ~ Easting, image = TRUE)

par(mfrow=c(2,3))
image(gam1,Northing ~ Easting)
image(gam2,Northing ~ Easting)
image(gam3,Northing ~ Easting)
image(gam4,Northing ~ Easting)
image(gam5,Northing ~ Easting)

See https://www.r-bloggers.com/r-as-gis-for-ecologists/ for plots https://www.r-bloggers.com/colored-3d-map/

References

Seemingly unrelated regression

Perperoglou, Aris, Willi Sauerbrei, Michal Abrahamowicz, and Matthias Schmid. 2019. “A Review of Spline Function Procedures in r.”